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Abstract 

We present a proposal for the full phase space distribution of the Milky Way halo. The model is 
axially and reflection symmetric and its time evolution is self-similar. It describes the halo as a set 
of discrete dark matter flows with stated densities and velocity vectors everywhere. We first discuss 
the general conditions under which the time evolution of a cold collisionless self-gravitating fluid is 
self-similar, and show that symmetry is not necessary for self-similarity. When spherical symmetry 
is imposed, the model is the same as described by Fillmore and Goldreich, and by Bertschinger, 
twenty-three years ago. The spherically symmetric model depends on one dimensionless parameter 
e and two dimensionful parameters. We set e = 0.3, a value consistent with the slope of the power 
spectrum of density perturbations on galactic scales. The dimensionful parameters are determined 
by the Galactic rotation velocity (220 km/s) at the position of the Sun and by the age of the Galaxy 
(13.7 Gyr). The properties of the outer caustics are derived in the spherically symmetric model. 
The structure of the inner halo depends on the angular momentum distribution of the dark matter 
particles. We assume that distribution to be axial and reflection symmetric, and dominated by 
net overall rotation. The inner caustics are rings whose radii are determined in terms of a single 
additional parameter jmax- We summarize the observational evidence in support of the model. 
The evidence is consistent with jmax = 0.18 in Concordance Cosmology, equivalent to jmax.oid = 
0.26 in Einstein - de Sitter cosmology. We give formulas to estimate the flow densities and velocity 
vectors anywhere in the Milky Way halo. The properties of the first forty flows at the location of 
the Earth are listed. 
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I. INTRODUCTION 



It has been established from a variety of observational inputs that approximately 23% 
of the energy density of the universe is in the form of "cold dark matter" (CDM)[l|. The 
CDM particles must be non-baryonic, coUisionless, and cold. "Cold" means that their 
primordial velocity dispersion is small enough that it can be set equal to zero for all practical 
purposes when discussing the formation of large scale structure. "CoUisionless" means that 
the particles have negligible interactions other than gravity. A central problem in dark 
matter studies is the question how CDM is distributed in the halos of galaxies, and in 
particular in the halo of our Milky Way galaxy. Indeed, knowledge of this distribution is 
essential for understanding galactic dynamics and for predicting signals in direct and indirect 
searches for dark matter on Earth. 

Galactic halos are coUisionless fluids and must therefore be described in 6-dimensional 
phase space. A full description gives the phase space distribution f{r,v; t) of the dark matter 
particles in the halo, i.e. their velocity (v) distribution at every position r. t is time. An 
important simplification occurs in the case of halos built of cold coUisionless dark matter 
because CDM particles lie in phase space on a thin 3-dimensional hypersurface. This fact 
implies that the velocity distribution is everywhere discrete [2j] and that there are surfaces 
in physical space, called caustics, where the density of dark matter is very large. See Fig. 
1 for an illustration. It is argued in ref. ^ that discrete flows and caustics in galactic 
halos are a generic and robust prediction of cold dark matter cosmology, even after all 
possible complications and reasons for skepticism have been considered. The reader may 
wish to consult ref. [3] for background information and a list of references. Discrete flows 
and caustics are seen in N-body simulations of large scale structure formation when special 
techniques are used to reach adequate resolution in the relevant regions of phase-space 

Bat 

The caustic ring model is a proposal for the full phase space distribution of the Milky 
Way halo. The model is axially symmetric, reflection symmetric, and self-similar in its 
time evolution. The model is an elaboration of the spherically symmetric self-similar model 
originally proposed by Fillmore and Goldreich [7], and by Bertschinger [sj. Our main purpose 
in presenting the caustic ring model is to enable its further comparison with observation. 
The observational evidence which we claim in support of the model thus far is listed in 
Section IV. 

Because the original self-similar halo model 0, [1] is spherically symmetric, one may be 
tempted to suppose that self-similarity and spherical symmetry are linked somehow, that 
the latter is required by the former. This is not so, however. In Section II, we discuss the 
general conditions under which the time evolution of a cold coUisionless self-gravitating fluid 
is self-similar. We show that self-similarity does not require symmetry. We consider three 
cases: no symmetry, spherical symmetry, and axial symmetry. In each instance, we derive 
the equations of motion for the reduced variables, and the initial conditions required by 
self-similarity. 

The spherically symmetric infall model is inadequate to describe the inner parts of a cold 
dark matter halo. Indeed, in a spherically symmetric halo, CDM particles necessarily have 
vanishing angular momentum. Each particle moves on a radial orbit and passes through 
the center at each throughfall. The center is then a caustic point where the density d{r) 
diverges approximately as ^, where r is galactocentric distance. In realistic halos, the dark 
matter particles carry angular momentum and therefore miss the center. There are still 
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inner caustics in that case [9j but they are spread on surfaces with cusps. The catastrophe 
structure of the inner caustics j9| depends on the angular momentum distribution of the 
infalhng CDM particles. If the angular momentum distribution is dominated by net overall 
rotation, the inner caustics are rings (closed tubes) whose cross-section is a section of the 
elliptic umbilic (-D-4) catastrophe [lo|, These objects are usually referred to as "caustic 



rings of dark matter" . Figs. and E] describe the cross-section of a caustic ring of dark 
matter and indicate what is meant by the ring radius a and by its transverse dimensions 
p and q. In the opposite case, where the velocity distribution is irrotational (V x -y = 0), 
the inner caustic has a tent-like structure described in detail in ref. [9]. Because there is 
evidence for caustic rings of dark matter (as opposed to tent-like caustics), in the Milky Way 
[l2| and in other isolated spiral galaxies [13], we assume that the angular momentum 



distribution is dominated by net overall rotation. 

The spherically symmetric self-similar infall model 0, @] was generalized in ref. 14] to 
include angular momentum of the infalling particles. The main purpose of that paper was to 
estimate the densities and speeds of the dark matter flows at a distance of 8.5 kpc from the 
Galactic center. It was shown in ref. |14] that the angular momentum distribution can be 
chosen so as to preserve self-similarity. For the sake of convenience, the model of ref. [3] was 
artificially made spherically symmetric by averaging an actual halo with angular momentum 
over all orientations. No such averaged halo can be physically realized in cold dark matter 
cosmology. However, the model of ref. [13] is adequate for its purpose of estimating the 
speeds and average densities of the dark matter flows at 8.5 kpc from the Galactic center. 



The model of ref. [ij] has inner caustics but they are spread over spheres. The model of 



ref. [ij] is an important precursor to the model presented here. The latter is essentially the 
former without the averaging over all orientations. 

It was assumed in ref. [IJ] that the universe is Einstein - de Sitter, i.e. spatially flat and 
matter dominated [VL = fim = 1). This is a convenient assumption when building a self- 
similar model because self-similarity requires the absence of a special time. The Einstein - 
de Sitter cosmology doesn't have a special time and is the only cosmology in which structure 
formation can be strictly self-similar. Subsequently, the accelerated expansion of the universe 
was discovered p^] and the Concordance Cosmology {Q^a = 0.27 and = 0.73) established 
[l|. The Concordance Cosmology is inconsistent with exact self- similarity since it has a 
special time, namely the time a few Gigayears ago when the universe passed from matter 
domination to dark energy domination. However, the Milky Way halo was assembled almost 
entirely before the onset of accelerated expansion, when the universe was well described 
hy Q = f^m = 1 cosmology. The breaking of self-similarity implied by the Concordance 
Cosmology has only a small effect on the formation of our Galaxy and will be ignored here. 

Nonetheless, the discovery of the accelerated expansion of the universe and the advent of 
Concordance Cosmolgy did affect the self-similar infall model in an important way because 
they led to a precise determination of the age of the universe. In ref. ij] it was assumed, 
consistently with Einstein - de Sitter cosmology, that the present age is related to the present 
Hubble expansion rate Hq by to = 3^ = 6.52 Gyr, where h is defined as usual hj Hq = h 
100 km/(s Mpc). Since the universe is now known to be 13.7 Gyr old [H, [iH], the results of 
ref. [13] remain valid provided only that h is replaced everywhere in that paper by ^'^^^^'^ 
= 0.476. 

The outline of our paper is as follows. In Section II, we discuss the general conditions 
under which the time evolution of a cold coUisionless self-gravitating fluid is self-similar. We 
show that spatial symmetry is not necessary for self-similarity. We discuss three cases - no 
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symmetry, spherical symmetry and axial symmetry - and obtain the equations of motion 
for the reduced variables in each case. In Section III, we present the model. We discuss 
its underlying assumptions, and describe its outer and inner caustics. We provide interpo- 
lating formulae to estimate the flow velocities and densities between the inner and outer 
caustics. We list the estimates for the inner and outer caustic properties which follow from 
the application of the method of adiabatic invariants to the model. In Section IV, we list 
the evidence in support of the model. In Section V, we list the densities and velocities of 
the first forty fiows on Earth. Section VI provides a summary. 



II. SELF-SIMILAR COLD DARK MATTER HALOS 

The purpose of this section is to obtain the conditions under which the time evolution 
of a cold self-gravitating coUisionless fluid is self-similar. We show that symmetry is not 
necessary for self-similarity. First, we discuss self- similar coUisionless fluids. Next, we discuss 
cold coUisionless fluids. Thirdly, we discuss cold and self-similar coUisionless fluids. In each 
instance we obtain the relevant equations, without assumptions of symmetry. Finally, we 
specialize to the cases of spherical and axial symmetry. 



A. General conditions for self-similarity 

Let /(r, v; t) be the phase space density of a coUisionless self-gravitating fluid at time t. 
It satisfles the coUisionless Boltzmann equation 



where 



$(f,t) = -Gm I d'r' I d'v ^[l'''": ^} (2.2) 



is the gravitational potential, m is the particle mass. 

Let us flrst discuss simple rescalings of the phase space distribution and evolution rate of 
such a system. Given that /(r, v; t) is a solution of Eqs. (12.11) and (12. 2p . consider 

nr,v-t) = AfC-,^-^-) , (2.3) 

where A, R, V and T are constants. It is straightforward to show that f'{r,v]t) is also a 
solution provided 

V = ^ and A = ^ . (2.4) 

Therefore any solution can always be rescaled in two ways: it can be made older (and hence 
more slowly evolving) and it can be made larger. The age t sets the scale for the density 
p since Gp ~ ^. The size R then sets the scale for the velocities f ~ y, and for the total 

mass M ~ pR'^ ~ Next let us consider self-similar evolution of the fluid, which is an 
invariance under time-dependent rescalings. 

The evolution is called self-similar if, in all aspects, the fluid remains identical to itself 
except for an overall rescaling of its phase space density by a factor A{t), of its size in the 
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spatial directions by a factor R{t) and of its size in the velocity directions by a factor V{t). 
Thus, the phase space density of a fluid with self-similar evolution satisfies the ansatz 



/(f.,7;«) = ^Wf^_._J (2.5) 

where F{x, ^) is a rescaled time-independent phase space density. The gravitational poten- 
tial is then 



W.t) = A{t)R{tYV{tY^[^-^^j (2.6) 

with 

^>{X) = -Gm [ d'x' [ d'u . (2.7) 

J J \x-x\ 

Substituting Eqs. (12. 5p and (12.61) into Eq. (12.11) . one finds that the self-similarity ansatz can 
be satisfied only if 

R{t)^t^ , V{t) = ^ and Mt) = ^ ■ (2.8) 
The rescaled phase space distribution must satisfy 

(1 - 3/3)F(x, u) + {V- m ■ II + ((1 - - f|) ■ 1^ = ■ (2-9) 

Thus, the evolution of a self-similar self-gravitating collisionless fluid is obtained by choosing 
l3 and solving simultaneously Eqs. (12. 9p and (12. 7p . with appropriate boundary conditions. 
Next, we consider a collisionless fluid which is self-gravitating and cold, but not necessarily 
self-similar. 



B. Phase space distribution of cold dark matter 

Cold collisionless dark matter (CDM) particles, such as axions or WIMPs, lie on a thin 3- 
dim. hypersurface in phase space. Indeed, at a sufficiently early time tin, all CDM particles at 
a given location r have the same velocity t''in(^, up to a small species-dependent primordial 
velocity dispersion 6v. We call the 3-dim. hypersurface the "phase space sheet". The 
thickness of the phase space sheet is the primordial velocity dispersion 6v. For CDM 6v 
is small, i.e. the sheet is thin [ll']. The number of particles is huge (approx. 10®'* axions 
and/or 10^^ WIMPs per galactic halo), so that the sheet is continuous. As time goes on, 
the sheet folds in phase space. 

To describe the evolution of the phase space sheet we label the particles by a continuous 
parameter a = (ai, ^2, as). The position of particle labeled a at time t is x(a, t). The parti- 
cle has velocity v{a, t) = ^{a, t). The phase space sheet is the set of points (x(a, t),v{d, t)) 
for all d. It has a time-independent mass density ^(5) in parameter space. At any point 
r in physical space, there is a discrete set of flows. The number A^(r, t) of flows at position r 
and time t is the number of solutions dj{f, t), j = 1, 2...N, of the equation r = x{d, t). The 
j**^ flow has velocity 

Ox 

^A^^t) = l^i^,ir,t),t) (2.10) 
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and density 

, , ^ , dM , ^ , ^ , , 1 , , 

<i,(r.t) = ^(a.(r.t))^^p^^^ (2.11) 

where 

D{a,t) = det(^^^^^ . (2.12) 

Caustics are located where D{dj{f,t),t) = for some j. A possible parametrization is 
a = Tin where r'in is particle position at an arbitrarily chosen but sufficiently early initial 

time tin. Other parametrizations may be more convenient, however. The flow velocities 
Vj{f,t) and densities dj{f,t) are parametrization independent. 
Since the particles are coUisionless, 

= -V[$(f,t) + $b(r,t)]|,=,K,) (2.13) 

where 

""''"d,(r'.i) _ ^ [ adM,^^ 1 



<t(r,t) = -G f Sr' Y: '^^=-Gf 



is the gravitational potential due to dark matter and $b is the contribution of ordinary mat- 



ter. The CDM phase space distribution is determined by specifying ^(a) and the positions 
rin(a) and velocities Vi-a^oi) of the particles at an initial time tin, and solving Eqs. fl2.13p and 
fl2.14p simultaneously. $^(7% t) has to be supplied by hand or derived in some way. 

We wish to apply the above to the construction of galactic halos. Galactic halos grow 
by accreting more and more of the surrounding dark matter. A convenient choice for the 
parameters a is as follows. Define a sphere of radius R{t) which grows with the halo. R{t) is 
chosen such that each particle traverses the sphere only once, in the inward direction. The 
particles are labeled by a = (6^, ip, r) where (6*, ip) are the spherical coordinates of the point 
where the particle crosses the sphere and r is the time when it does so. For given -R(t), the 
model halo is specified by giving the velocities 'yin(6', r) of the particles when they cross 
the sphere and their infall rate per unit solid angle: 



dVtdr sin 9 dOdipdr 

Note that the number of degrees of freedom in ViniO, ip, r) is slightly reduced by the freedom 
to choose R{t). For example, we may choose R{t) such that the particle which is at the 
intersection of the sphere with the x-axis is at ffist turnaround, for all t. In that case, 
'Uin(f , 0, r) ■ X = for all r. 

The gravitational potential due to dark matter may be separated into two parts: $ = 
<l>in + $cx, with 

"* dM 1 



M = -Glml (2.16) 







the contribution of particles that have already fallen through the sphere, and 

_ dM 1 
dVtdr |r — x{6, (p, r; t)\ 



°° dM 1 
$^^(f,t) = -G I dQ I dr^^——^ (2.17) 
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the contribution of particles that have not yet done so. In the spherically symmetric case $ex 
can be ignored since it is r-independent inside the sphere. In the non-spherically symmetric 
case, to obtain $ex, the equations of motion have to be - in principle - integrated backward 
in time to determine the positions of the particles outside the sphere. This procedure would 
then include all the particles in the universe. In practice, the contribution of particles outside 
the sphere may be given by hand, or ignored. 



C. Self-similar evolution of cold dark matter 

In this subsection, we set = because the presence of ordinary matter generally spoils 
exact self-similarity by introducing special time and special length scales. However later, 
when discussing our Milky Way halo model, we will add the baryon contribution to the 
gravitational potential by hand, maintaining self-similarity, as described in Section III. 

We saw in subsection II. A that a self-similar halo grows in physical space as a power law: 
R{t) oc t^. Hence we choose the sphere introduced in subsection II. B to have radius 



R{t) = Ro\^-j (2.18) 

where to is an arbitrary reference time. Under the self-similarity assumption, the particle 
that falls in at time r and location {6, (p) moves in exactly the same way as the particle that 
falls in at (6*, ip) at any other time r', except that its trajectory is rescaled in physical space 
by and in time by ^. Hence there exists a function \{9, p)\ ^) such that 

x{e,p>,T-t) = R{T)\{e,^-i = -) . (2.19) 

T 

Also the infall rate ^^(^, f, t) must have a specific time dependence. The total dark matter 
mass within the sphere of radius R{t) is 

Mit) = f dVL f dT^^{e,ip,T) =m [ d^r [ d^vf{r,v;t) 
J Jo dildr Jr<R{t) J 

'IxKi 

where we used Eqs. (12. 8p . Hence we need 



mA{t)^^ f d\ I d''uF{x,V)(xe^-^ , (2.20) 



dM ,^ , /^\3(/3-l) 



dQ.dT \to J dVtdr 

Finally, the initial velocities must be of the form 

R{t 



(e,^) . (2.21) 





yin{e,^) (2.22) 



to be consistent with Eq. (I2.19p . 

Substituting Eqs. flCTl) and ([221]) into Eq. (EUD, we obtain 



(2.23) 
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with 



,^0/ Jo ^3^-1 J dndT 
With $b = 0, the equations of motion (I2.13P become 



(9,^)-^^^ . (2.24) 



The initial conditions are 

{e,^-i) = u,^{e,^) , (2.26) 

where h{6,(p) is the unit vector in the direction {6,ip). The halo is determined in terms of 
P, Vin{0,v) and J^|o(^,<^) by solving Eqs. ^I?2Mi and fl2:25|) with initial conditions ^Mj- 
Again, the freedom to choose R removes one degree of freedom from i?in(6', v?)- For example, 
R may be chosen such that z7in(|, 0) ■ x = 0. 



D. The spherically symmetric case 

For a halo of cold collisionless particles to be spherically symmetric, the particle orbits 
must be radial, i.e. each particle has zero angular momentum. Thus 

z7i,(^, ^) = Ui^h{9, ^) . (2.27) 

Furthermore z/in and -^^lo must be independent of {O^ip). We choose R{t) to be the 
turnaround radius at time t, so that = 0. We have then 



Because of spherical symmetry, f^(x) gsts contributions only from ^ such that |^A(^) < x- 
In particular, for x < 1 there is no contribution to ^ 
may write therefore 



In particular, for x < 1 there is no contribution to 4^ from the integral over < ^ < 1. We 



-T^ix = X^ 2.29 
OX 8 X 



where 

3 



The factor ^ is introduced so that M{1) = 1 [see Eqs. dOHil and (IQIjl below]. The 
equation of motion is 

M[x = -^] ■ (2.31) 
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Equations fl2.30p and fl2.3ip were first obtained by Fillmore and Goldreich |7|. In terms of 
the parameter e which these authors introduced, = | + ^- Bertschinger [81] discussed the 
case e = 1 (/? = §)• 

Following these authors, consider a spherically symmetric overdensity in an otherwise 
homogeneous, fiat and purely matter dominated universe. Assume the overdensity has a 
power law profile 

^ « hrr 2.32 

Mi \mJ ^ ^ 

where Mj is the mass interior to radius rj in the homogeneous universe, and SMi is the 
corresponding overdensity. The exponent e must be in the range 0<e<l. e = l describes 
the case where the overdensity is a point mass. It can be shown 0, U 14| that the time 

2 2 

evolution of such an overdensity is self-similar with scale factor R(t) oc ta + g?. The mass 
interior to the turnaround radius R{t) is 

2 _ 

Note that M{t) oc ts?, consistent with Eq. (12.201) . Therefore, Eq. (12.301) may be rewritten 

MX) 4 r -^e (, - ■ (2.34) 

We have = 1 since the argument of the B-function in Eq. (I2.34p is positive for all ^ 

when X = 1- Thus 

M{r,t) = M{t)M ( ) (2.35) 



m 



is the mass interior to radius r at time t. 



E. The axially symmetric case 

In the axially symmetric case, the model is specified by giving /3, ^^|o('^) and the initial 
velocity field 

ViniO, (fi) = f Uria) + 6 ueia) + (p v^{a) , (2.36) 

where a = | — 0. If reflection symmetry is assumed in addition to axial symmetry, ^^|o) 
Vr and are even under a —a, whereas is odd. One degree of freedom is removed 
from z7in by the choice of R. For example, we may choose R such that all the particles in 
the equatorial plane are at first turnaround, in which case //^(O) = 0. 
In view of the symmetry, we have 

K^^V'^i) = zX-Xa,i) + \p{a,i)[xcos-i{a,ip]i) +ysin-i{a,ip]i)] , (2.37) 

where the (p dependence of 7 is trivial, i.e. 

7(a,<^;0 = <^ + A(a,0 • (2.38) 
The z component of angular momentum is conserved: 

j^{a) = \p{a,if — {a,C) = v^{a) cos{a) . (2.39) 
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A{a,^) is determined by solving Eq. (12.390 with the initial condition A(a,0) = 0. The 
equations of motion for the other degrees of freedom are 



d^X, _ _ 1 9^ 



X 



with initial conditions 



\z{a,l) = sina , Ap(a, l) = cosa 



(a, 1) = sin a Vr\(y) — cos a ve^a) 
(a, 1) = cosa z/j,(a) + sina z/e(Q;) . (2-41) 



III. THE MODEL 

A. Assumptions and general properties 

To first approximation when discussing its large scale properties, the model is spheri- 
cally symmetric. So we start with the spherically symmetric self-similar infall model 0, [1] 
introduced in Section II. D. The dimensionless properties of that model are determined in 
terms of a single parameter e = -^^zq- cin average sense, e is related to the slope of the 



power spectrum of matter density perturbations [l6j. On galactic scales (Mpc), e is of order 



0.3 14 



In the spherically symmetric model, the angular momentum of each particle vanishes and 
all orbits are radial. The particle positions are given by 

r{e,^,T-t) = R{T)h{e,^)x(^-^ , (3.1) 

where R{t) is the turnaround radius at time t, and r is the time when the particle reached 
first turnaround. The function A(^ = ^) is obtained by solving Eqs. (12.301) and (I2.3ip . This 
can be done numerically on even a modest computer. Fig. [2] shows the function A(^) for 
e = 0.3. Fig. [3] shows the corresponding phase space distribution. 

To apply the description to an actual halo, we must provide two dimensionful properties 
of the halo. Two of the best measured properties of the Milky Way halo are the rotation 
velocity at the position of the Sun Vrot^fQ = 8.5 kpc) = 220 km/s, and the age to = 13.7 
Gyr. We will ignore the uncertainties on these two numbers. 

Fig. m shows the dimensionless rotation velocity squared z^^(e, x) a function of dimen- 
sionless galactocentric distance x = for several values of e. z/(e, x) is defined by 



2 , , GM(r) 2. r^GM . 
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where M(r) is the mass interior to r, and M = M{R). Therefore 



Mr X 

M{e,x) is given by Eq. fl^TMD . Note that z/(e, 1) = 1. Combining Eqs. ^M) and ([S2D, we 
have 



R = —v,otir) t ^) . (3.4) 

71 R 

The turnaround radius R is obtained by solving Eqs. (13. 4p with t = 13.7 Gyr and frot('"0) 
= 220 km/s. Table [I gives R and M for various values of e. This is an update of Tabled in 
ref. 0. 

The model is characterized by a set of outer caustics located where the outflows reach 
their maximum radii Rn {n = 1,2,3...) before falling back in. The properties of the outer 
caustics will be discussed in detail in Section III.B. In addition to these outer caustics, a 
galactic halo necessarily has inner caustics [9|]. In the spherically symmetric model, the inner 
caustics have all collapsed to the central point. Indeed, since its angular momentum vanishes, 
each particle passes through the center in the course of each in and out fall. The center 
is then a singular point, where the density diverges. Fig. H] shows that, in the spherically 
symmetric model, the halo contribution to the rotation curve is constant near r = for e ~ 
0.3, implying that the density d{r) oc ^ as r 0. 

The spherically symmetric model does not adequately describe the inner parts of the halo. 
In actual galaxies, the central parts are dominated by baryons and the halo contribution to 



the rotation curve goes to zero as r — >■ 0. As pointed out in ref. [l4|, the depletion of the 
inner halo, as compared to the prediction of the spherically symmetric self-similar model, is 
neatly accounted for by the angular momentum that the dark matter particles are expected 
to carry because angular momentum causes the dark matter particles to avoid the center. 
Indeed, by constraining the distance scale (effective core radius) over which the inner halo 
is depleted, one can estimate the average specific angular momentum of the halo particles 

Q. 

When angular momentum is introduced, the inner caustics spread over surfaces with 
cusps. (The word "cusp" has in this paper the same meaning as in Catastrophe Theory, 
unrelated to the notion of "cuspy halos".) In ref. the catastrophe structure of the inner 
caustics was investigated as a function of the angular momentum distribution of the infalling 
particles. If the angular momentum distribution is dominated by net overall rotation, the 
inner caustic is a ring (i.e. a closed tube) whose cross-section is a section of the elliptic 
umbilic (D-a) catastrophe jll]. We call this type of inner caustic a "caustic ring of dark 
matter" [lol.]ll|. Figs. [5] and E] show such a ring in cross-section. Fig. [5] shows the flow of 
dark matter in the neighborhood of the caustic, whereas Fig. [6] shows the definition of its 
radius a and transverse sizes p and q. Caustic rings of dark matter are described in more 
detail in Section III.C below. When the initial velocity distribution is dominated by a curl- 
free (non rotational) component, the inner caustic has a "tent-like" structure, which may 
be described qualitatively as a caustic ring of dark matter whose (inner) radius a shrunk to 
zero while its outer radius a + p was held fixed. See ref. [9| for a detailed description. 

As summarized in Section IV, there is evidence for caustic rings of dark matter in the 
Milky Way and other spiral galaxies. The evidence is for caustic rings, as opposed to the 
tent-like inner caustics that occur when the angular momentum distribution is irrotational. 
Hence we assume that the angular momentum distribution of the infalling dark matter is 
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dominated by net overall rotation. We assume furthermore that it is axially and reflec- 
tion symmetric, and consistent with self-similarity. To our previous spherically symmetric 
description depending on one dimensionless parameter e, we add therefore the dimension- 
less initial velocity field of Eq. fl2.36p . The specific angular momentum of the particle at 
turnaround radius -R(r) and angular coordinates (6*, ip) is then 

r(9,^,T) = ^ji9,^) (3.5) 
r 

where 

— * ^ 

j{e,(f) = h{9,ip) X Ui^{9,ip) = -u^{a)9 + ug{a)ip . (3.6) 

The assumption that the initial velocity distribution is dominated by net overall rotation 
means z/^^ >> |z/r|, |z^6»|- The simple case of initial rigid rotation is ^^{a) = z/<^(0)cosa, 

= 1^0 = 0. 

It was mentioned already that, when angular momentum is added to the spherically 
symmetric infall model, the inner parts of the halo are depleted and the halo contribution 
to the rotation curve goes to zero when r 0, consistent with observation. In the Milky 
Way and in many other large spiral galaxies, the baryon contribution to the rotation curve 
is such as to make the curve approximately fiat at small r, and with approximately the same 
value at small r as at large r. In other words, the baryons and the dark matter "conspire" to 
keep the rotation velocity approximately constant for all radii where it has been measured. 
Now, recall that the rotation curve of the spherically symmetric self-similar model is very 
nearly constant for all r when e ~ 0.3; see Fig. HI Following ref.[l3|, we take account 
of the gravitational force exerted by the baryons in the Galaxy by using the potential 
of the spherically symmetric model to calculate the motion of the dark matter particles 
in the axially symmetric model. This should be a good approximation when calculating 
the properties of low n flows since they are made of particles falling in from distances 
much greater than the size of the baryonic disk. Furthermore, the caustic rings are made 
of particles that move in or near the equatorial plane and therefore respond to the same 
gravitational forces as produce the nearly fiat rotation curve of the Galaxy. 

In ref. the average value j of j{9, ip) = \j{9, Lp)\ was estimated to be 0.2 by requiring 
the halo contribution to the rotation curve to approach zero as r — > over a distance scale 
(effective core radius) consistent with observations. Subsequently, evidence was found for 



caustic rings of dark matter in the Milky Way |10l . |12|. This evidence is consistent with 
the maximum value jmax of j{9, ip) being equal to 0.263. jmax is to be identified with z/<^(0). 
Although the precise relationship between j and jmax depends on the full j(6', (y9)-distribution, 
the values j ~ 0.2 and jmax — 0.263 are consistent with one another. For example, for the 
simple case of initial rigid rotation, one has jmax = and hence jmax = 0.255 for j = 0.2. 

Both jmax — 0.263 and j ~ 0.2 are the outcome of fits to observations in the old fl = Q^a = 
1 cosmology with Hq = 70 km/ (s Mpc) in which the age of our galaxy is 9.31 Gyr. Since we 
now adopt the age of 13.7 Gyr, while keeping the rotation velocity (220 km/s) unchanged, all 
length scales are stretched by the factor 13.7/9.31= 1.47. Because the effective core radius 
is approximately proportional to j and the caustic ring radii approximately proportional 
to jmax, the values of these parameters should be reduced by the factor 1.47 compared to 
the old fl = flra = 1 model. So the fit to observations yields j ~ 0.14 and jmax — 0.179 
in Concordance Cosmology. More precisely, because the caustic ring radii are not exactly 
proportional to jmax, the fit of the model to the evidence in Concordance Cosmology (see 
Section IV.D) yields jmax = 0.186, which is the value adopted in this paper. To avoid 
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confusion, we write henceforth jmax,oid for values in the old Q = = 1 model, and simply 
irnax for values in the Concordance Cosmology. 

For the sake of definiteness, the model presented here assumes e = 0.3. 

B. Outer caustic spheres 

Caustics are at the boundaries between regions in physical space with differing number of 
flows. On one side of a caustic there are two more flows than on the other. When an isolated 
halo is approached from the outside, the local number of flows increases. First, there is one 
flow, then three, then flve, and so on. The boundary between the region with one flow and 
the region with three flows is the location of the first outer caustic. Likewise, the second 
outer caustic is at the boundary between the region with three fiows and the region with 
five flows, and so on. Sec Fig. 1 for an illustration. 

The outer caustics are simple fold (742) castastrophes located on topological 2-spheres 
surrounding the galaxy. The density in the 2 extra flows on the inside of an outer caustic 
diverges on the 2-sphere, with the following characteristic behaviour: 



where x is the distance to the 2-sphere, measured positively on the inside, Q{x) is the step 
function, and A is a constant which we call the fold coefficient. In the absence of symmetry, 
A varies with location on the 2-sphere. The divergence of as x — > 0+ is cut off by the 
velocity dispersion of the dark matter particles. 

Let us emphasize that the existence of outer caustics with the properties listed in the 
previous paragraph follows exclusively from the existence of cold collisionless dark matter. 

1. Radii Rn 

Because angular momentum has little influence on the outer caustics, we use the spheri- 
cally symmetric model to describe their properties and, in particular, to estimate the galac- 
tocentric radii i?„ where they occur and their fold coefficients An- As before, each spherical 
shell is labeled by the time r when it was at its first turnaround. At time t, the radius of a 
shell is 




(3.7) 




(3.8) 



where we used R{t) oc t^. The outer caustics are at radii 




(3.9) 



where the ^„ are the locations of the maxima of A(^) = ^ ''A(^). 
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2. Fold coefficients 



From Eqs. fl2.34p and fl2.35p . one may derive the following formula for the density d{r,t) 
as a sum over flows 



1 M{t) 
67re i?(t)r2 



1 



(3.10) 



where the ^j(r, t) are the solutions of r = R(t) A{C,). Near a maximum ^„ of A(^), we have 



(3.11) 



with ^^{^n) < 0. Let us label j = {n+) and (n—) the two flows which form the caustic at 



r = Rn- Then (the t dependence is suppressed in the remainder of this subsection) 



J] 



(3.12) 



for r close to, but less than, Combining Eqs. (13.101) and (13.121) . we find 

1 An 



with 



dn±ir) 



1 M 



QiRn~r)[l + Oi^/Rn-r)] 



(3.13) 



(3.14) 



The 0{y/Rn — r) remainder in Eq. (13.13^ has opposite sign for (i„+(r) and dn-{r) so that 

A, 



dn{r) = dn+{r) + t/„_(r) = ^^=e(i?„ - r)[l + 0{Rn - r 



(3.15) 



Table [TTl lists the present Rn and An values in the spherically symmetric infall model with e 
= 0.3, for n = l,2, ... 20. 



C. Inner caustic rings 

It was mentioned already that dark matter halos necessarily have inner caustics, that the 
catastrophe structure of the inner caustics depends on the angular momentum distribution 
of the infalling dark matter particles [9|], and that the inner caustics are rings if, as we assume 
to be the case, the angular momentum distribution is dominated by net overall rotation. 



The caustic ring singularity was described in ref. [llj which the reader may wish to consult 
for background information. 
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In the limit of axial and reflection symmetry and where the transverse sizes p and g of a 
caustic ring are much smaller than its radius a, the distribution of dark matter particles in 
the vicinity of the caustic ring is given by the particle positions 

z{a,ri) = bar] 

p{a,r]) = a+^u{r] -r]of -^sa'^ . (3.16) 

We use cylindrical coordinates {z, p, ip) for position in physical space. Eqs. (13.161) give particle 
positions at a particular time t, which is not shown explicitly. The particles are labeled by 
parameters {a, rf). As before, a = ^ — 9 where 9 is the polar angle of the particle at the time 
of its first turnaround, rj is the time when the particle crosses the z = Q plane in the course 
of its flow through the caustic ring, t — rj can be thought of as the age of the particle. The 
particles labeled (a,?]) form a circle of radius p{a,ri) at a height z{a,r]) above the z = 
plane. 

Fig. Oplots (p(a, T]), z{a, t])) for continuous rj, and discrete values of a. The lines in Fig.O 
are the trajectories of the particles forming the flow, except that positions are plotted as a 
function of age, whereas for ordinary trajectories position is plotted as a function of time. 
Let us call the lines of Fig. [5] "age trajectories". Fig. [5] shows that particle density diverges 
on a closed line which has the shape of a isosceles triangle, but with cusps instead of angles. 
We call that shape a"tricusp". The location of the tricusp is shown in Fig. [6] for the flow of 
Fig. [5l It is the envelope of the age trajectories. There are four flows everywhere inside the 
tricusp and two flows everywhere outside. The caustic, i.e. the surface where the density 
diverges, lies at the boundary between the region with four flows and the region with two 
flows. 

The physical space density is given by [ll| 



,, , 1 dM cos a . . 

d{p,z) = - } ipri- («,'?) 77^7 77 (3-17) 



dVtdrj ' \D{a,'r])\ 



where aj{p,z) and rij{p,z), with j = I ... N{p,z), are the solutions of p^oi^rf) = p and 
z{a, Tj) = z. N{p, z) is the number of flows at position (p, z); thus, N = 4 inside the tricusp 
and N = 2 outside. £*(«, rj) is the Jacobian determinant of the map (a, rf) (p, z): 

D{a, ri) = det ( = -b[ur]{r] - r]o) + sa^] . (3.18) 

dM _ dM 



d{a,r]) 

— n J J is the mass faUing in per unit solid angle and unit time. The tricusp 

a\ldri 2tt cos adadrj o jr o r 

perimeter is the locus of points (p(a, r]), z{a, rj)) for which D{a, rj) = 0. We call p and q the 
sizes of the tricusp in the p and z directions respectively; see Fig. [61 They are given by 



1 2 v27 b 

2 A Wus 



Let us write the velocity of particle labeled (a, rj) as v{a, rj) = v^{a, ri)(p+Vp{a, ri)p+Vz{a, ri)z. 
The main component of velocity is in the if direction: ~ t>. In the case of a stationary 
flow, the velocity components in the p and z directions are 

dp dz 
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Here we use the fact that, in case of stationary flow, the particle positions are functions 
only of their age t — rj. Caustic rings grow in mass and radius on cosmological time scales. 
Therefore stationarity is not an exact property of caustic rings. However, it is an excellent 
approximation. In the self-similar model the growth in physical size is tantamount to ex- 
pansion in all directions by the scale factor R{t). Hence Vp in Eq. (13.201) should be replaced 

by + and Vz by Vz + ^z. However the corrections are small, of order 1% for n = 1 and 
less for n > 2, and will be ignored henceforth. For stationary flow, the speed v is related to 
u by 



u=— (3.21) 

a 

since u is the centrifugal acceleration of the particles at {z, p) = (0, a). 

Our description of a caustic ring in the limit of axial and reflection symmetry, and where 
the transverse sizes p and q of the ring are much smaller than its radius a, is in terms of six 
parameters: a, b, u, rjo, s and We now turn to the self-similar infall model to obtain 
predictions for many of these parameters. As was discussed in Section 111. A, the particles 
are assumed to fall in the gravitational potential produced by the mass distribution M(r, t) 
of the spherically symmetric self-similar infall model. The equation of motion is obtained 
therefore by combining Eqs. (I2.25P and (12.291) : 



de 

and using for A4{x) the solution of Eqs. (I2.30p and (12.311) . The initial conditions are 

dX 

x{e,v;i) = h{e,^) , —{e,^;i) = uue,v) (3.23) 



where z7in(6', ip) is the initial velocity defined in Eq. (12.361) 



1. Radii 

The radius of the n^^ caustic ring is the distance of closest approach to the galactic 
center of the particles in the equatorial plane for the rfi^ in and out flow. To obtain the 
an we solve Eq. (I3.22p with a = 6' — | = 0. By reflection symmetry, z^6»(0) = 0. Also, 
^'r-(O) = by definition of the turnaround radius R. Eq. (13.22^ implies conservation of 
angular momentum. We have 

A^(0 ^ = ^.(0) (3.24) 

for the particles in the equatorial plane. The quantity i^<^(0) was called jmax in previous 
publications. For the sake of consistency, we continue to use this name here. The equation 
of motion for the radial degree of freedom is then 

rfe " A3 8 ^ A2-^ ; ■ ^ ^ ^ 

The radial coordinate of the particles that reached first turnaround at time r is 

r(r,t) = R{t)\ (^jVnax, ^) = /?(t)^A(j^ax,0|^^, (3.26) 
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where A(jniax,0 is the solution of Eq. (13.251) with initial conditions: A(jmax, 1) = 
1) ^(jmax, 1) = 0. The caustic ring radii a„ are 

a„ = i?(t)A(j^ax,a (3.27) 

where the ^.^ are the locations of the minima of A(jmax,^) = 'C~^A(jmax, 0- The second 
column of Table UTTl lists the a„ for n = 1, 2, ... 20, in the e = 0.3 model with jmax = 0.186. 
This value of jmax in the Concordance Cosmology was determined from a fit of the a„ to a 
set of rises in the Milky Way rotation curve, as discussed in Section IV. 



2. Accelerations Un 

The particles at {z, p) = (0, a„) are moving in the direction with speed Vn = ^ where 

R^{t^) R^jt) 1 

Jmax Jmax T / 2/3—1 ' (3.28) 

is their specific angular momentum. Combining Eqs. (13.271) and (I3.28p . we obtain 



The Un parameters are then 



Un = ^ . (3.30) 

an 

The speeds f „ are listed in the third column of Table UTTl for e = 0.3 and jmax = 0.186. Note 
that the f„ only have a weak dependence on j^ax because the denominator in the RHS of 
Eq. (I3.29P is proportional to jmax in the small jmax limit. 



3. Infall rates ^ 



The infall rates are obtained by noting that M{t) oc ts^ and assuming that the infall is 
isotropic. We have 



dM 



dM .dr. 



dfldr] n dVtdr dr] 



1 2 M{t„) Tn 1 M(t) 1 



Atx 3e Tn t 



Gire t c'~e 



(3.31) 



f 3. 



That ^ 

drj 



dr I 
dt 



-- follows from self-similarity since, in the absence of any time scale, 
Tn and t must be proportional to one another. The fourth column of Table IIIII shows the 



dM 
dQdrj 



values in the model. 



17 



4- ??0,ni Sn 



The remaining parameters 6„, r/o.n and s„ are also predicted by the self-similar infall 



model but are related to more detailed properties of near a = 0. We have [11 



6„ = t;„[cos(5„(0) + 0'(O) sin5„(0)] (3.32) 

and 

_ an_ 0'(O) cos(5n(0) - sin(5„(0) 

0'(O)sin5„(O) + cos5„(O) ' ^ ' ' 

where 

0'(O) = --^^(O) (3.34) 



and 



Sn{0) = -^+Jm..J^ . (3.35) 



Ref. [llj gives also a formula for the s„ involving higher derivatives of 2/51(0;) and z/(^(a) at 
a = 0. Generally, s„ and a„ are of the same order of magnitude. 



D. Interpolating between inner and outer caustics 



The model predicts the full phase space distribution of the halo at all times, i.e. it predicts 
the number N{f, t) of flows at all locations, the flow densities dj{f, t) and the flow velocities 
Vj{f,t). To obtain these quantities one uses, in principle, the following procedure. Obtain 
\{9,ip;^) by solving Eqs. fl3.22p with initial conditions (I3.23P for all {9,ip). Then find the 
solutions {9,(p,T)j, j = 1,2...N, of 



R{t) [^yXie,^-/-] 



(3.36) 



The number of solutions is the number of flows at location r at time t. The velocities of 
the flows are 



Their densities are 



dj{r,t) 



R{t) (Tj\P-^ dX 



t 



dM 
dVtdr ' 



t 



t 



sin 6 



hfj;^ = — ) 



r 



, , / dix,y,z) 

' d{e,ip,T) 



(3.37) 



(3.38) 



It is certainly possible to carry out this procedure in practice, but it would be laborious. We 
note that the underlying formalism keeps track of all particles at all times. Since the particles 
are indistinguishable, this is far more information than the full phase space distribution at all 
times. Many practical tests of the model are sensitive only to the phase space distribution. 

The purpose of this section is to give approximate formulas for the full phase space 
distribution today, based upon the smooth interpolation between inner and outer caustics. 
The generalization to other times is straightforward since all lengths scale as t^, all velocities 
as t^~^ and all densities as t~^. 
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At an arbitrary location f which is not within the tricusp tube of a caustic ring, the 
number of flows is N{r) = 1 + 2 0(-Rn — f) where B(x) is the step-function. For example 
at r = 130 kpc, there are 9 flows because four outer caustics have radius larger than 130 
kpc; see Table [TTl If r is within the tricusp tube of a caustic ring, the number of flows is 
A^(r) = 3 + 2 6(-R„ — r) because there are two more flows inside the tricusp tube of a 
caustic ring than outside. 

For all r the flow speeds f„±(r) can be read off from Fig. [31 This neglects the effect of 
angular momentum on flow speed, a small effect. The flow speeds can also be estimated 
more crudely by using the formula 



Vn±{r)=v,ot\l2ln{^] . (3.39) 



r 



Eq. (13.39P would follow from energy conservation in the static gravitational potential $(r) = 
f^Q^ln(i). Because the actual potential is time- dependent, energy conservation is only 
approximate. It is least reliable for n = 1 since the gravitational potential is changing 
appreciably on the time scale of the first infall. 

To estimate the density of the n^^ in and out flows one may use 

provided r >> a„. Eq. fl3.40p expresses particle number conservation in the spherically 
symmetric case. The infall rates are given by 



dM 



dVtdt 



(,) = (3.41) 



where the in±ij') are the solutions of r = R{t)K{^). Eqs. fl3.4ip is obtained through the 
same steps as Eq. 03.311) . The infall rates are smooth functions of ^ and may be obtained 
by interpolating between their values at the inner caustics, given in Table IIIII 
In the spherically symmetric case, the flow equations near radius r are 

z{a,vi) = v{r) rj sin a , p(a,?7)=f(r) \rj\ cos a . (3.42) 

where v{r) is the speed of the flow at radius r. As before, rj is the time when the particle 
crosses the 2 = plane. These equations are valid for r » a where a is the radius of 
the inner caustic made by that flow. Close to the inner caustic, the flow is described by 
Eqs. f l3.16l) . Comparing Eqs. (13.161) and (I3.42p suggests an interpolating formula for the case 

h = V and s = a: 

z = V rj sm a , p= \/ ^ v^{rj — rj^y cos a . (3.43) 

2 

Indeed, since u = Eqs. (13.431) reduce to Eqs. (I3.16P for small rj and small a when b = v 
and s = a, and reduce to Eqs. (I3.42p for large rj and all a. To accomodate b ^ v and/or 
s a, one may generalize Eqs. (I3.43P to 

z = V rj sm{gz{rj)a) , p = a/o^ + f 2(77 - r/o)^ cos{gp{rj)a) (3.44) 

where gziv) 9p{v) smooth functions which approach respectively ^ and when 
?7 — > 0, but approach one when rj » rjo. 
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Sufficiently far from a caustic ring, one may choose to neglect its transverse dimensions 
p and q. This is done by setting r^o = 0. In that limit, Eqs. (13.431) may be easily inverted to 
obtain ri{z,p) and a{z,p). One can then express the density and velocity fields directly in 
terms of position : 



— Y - — + a/ (r^ — a^)^ + 4a^z^ 



=Fsign(2;)^-^y ^ (r"^ - + ^/(r^ - a^)^ + 4a2z2j + - ^^(H - 0^)2 + 402^2^ 



1 dM 1 , , 

d+ = d_ = --— = , 3.45 

V dndr ^(r2 - a2)2 + 4a2^2 

where the =F signs are for the down and up flows. 



E. Adiabatic approximation 



Except during the ffist couple throughfalls, the gravitational potential seen by a dark 
matter particle is slowly varying on the oscillation time scale of its radial coordinate. By 
exploiting the adiabatic invariant for the radial motion of the particle, it is possible to 
estimate its trajectory without resorting to numerical integration. This method was used 
in ref. [3] to show that the halo density c? oc ^ as r ^ in the spherical infall model with 



e < 



|. It was further used in ref. 



141 ] to obtain estimates of the densities and speeds of 



the flows at our distance from the galactic center. The adiabatic method gives predictions 
for the quantities of interest here. We state these predictions without proof, leaving their 
derivation to a future publication. 

For the outer caustics, the adiabatic method predicts 



Rn 
■An, 



4 4 + 3e 
9e 



n + 1 



4+3£ 



M 



37r2e R 



r6 



9e 
4+3e 



(3.46) 



To obtain these results, one assumes that the adiabatic invariant for radial motion is constant 
and that the potential is exactly logarithmic, i.e. $(r, t) = — frot(^)^ ln(r). 
For the inner caustics, the adiabatic method predicts 



Sn 



4 4 + 3e 

m 

7 

m 



[n 



4+3e 
9e 



t 



7 Jn 



9e 
4+3e 



(3.47) 
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where 7 is given in terms of jmax by 



To obtain Eqs. (13.471) . one makes the same assumptions as for Eqs. fl3.46p plus the assumption 
that jmax is smalL 

Comparison with the values listed in Tables [IT] and IIIII shows that the adiabatic method 
agrees with the results of numerical integration usually to within 10% or 20%. The overall 
agreement may be surprising considering that the first couple of throughfalls are not in the 
adiabatic regime, and that an exactly flat rotation curve is assumed to derive the predictions 
of the adiabatic approximation. Fig. H] shows that the rotation curve of the self-similar infall 
model is only approximately flat, even for e ~ 0.3. 

The adiabatic method is useful because it provides a check on the results of numerical 
integration. It explains patterns in the numerical results that are otherwise mysterious. In 
particular, it predicts that f„ is n-independent and that i?„ and a„ cx: ^ for large n; see 
Eqs. f l3.46l) and (13.471) . These properties are descriptve of the values listed in Table [ITl and 
IIIII Finally, Eqs. (13.461) and (13.471) enable quick estimates of all model properties for e 7^ 0.3. 



IV. OBSERVATIONAL EVIDENCE 

In this section, we list a set of observations which are consistent with the caustic ring 
halo model, and thus lend it support. 



A. Flat rotation curves 

The model predicts the Galactic rotation curve at large radii r all the way up to the 
turnaround radius R. The rotation curves for various values of e are shown in Fig. |3| 
in dimensionless units. They are nearly fiat when e ~ 0.3, implying that the density d 
is approximately proportional to ^ at large r all the way up to the turnaround radius 
R, whereabouts d reaches the average value of the cosmological dark matter density. In 



contrast, the density profiles predicted by computer simulations [17| behave as d{r) oc ^ at 
large r, implying that the rotation velocity decreases as 

The rotation cuve of the Milky Way is consistent with being flat up to the largest radii, 
of order 20 kpc, where it has been measured. More importantly, the rotation curves of spiral 
galaxies in general are flat at large radii up to the largest radii, of order 70 kpc, where they 
have been measured. A collection of extended well measured rotation curves was published 



in refs. UJ, 19 



Furthermore, studies of the dynamics of satellite galaxies i2(f\ provide evidence that the 
d oc behavior extends to r ~ 200 kpc. Weak lensing studies of the distortion of galaxy 
shapes by the gravitational fields of foreground galaxies ^] provide evidence that the be- 
havior extends to r ~ 370 kpc. 
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B. Effective core radii 



At small r, the Milky Way rotation curve is well accounted for by the baryonic matter 
in the Galactic bulge and disk. So, the halo contribution to the rotation curve should be 
suppressed at small radii and should vanish at r = [23] . The self-similar model generalized 
to include angular momentum of the dark matter particles [3] does have that property. In 
ref. [13], the 'effective core radius' rc^es of a galactic halo was defined as the radius at which 
the halo contribution to the rotation velocity squared is half its value, v'^^^, at large r where 
the halo dominates. rc,eff is proportional to the average amount of dimensionless angular 
momentum j. The value of r^efr for the Milky Way, estimated by modeling the bulge and 



disk contributions to rotation curve, implies jow ~ 0.2 [14 



C. Combined rotation curve 

Table UTTl gives the model predictions for the caustic ring radii a„. The predictions may 
be approximated by 

40 kpC / Wrot \ /jmax,old' 



, . , , , (4.1) 

n \220 km/s J \ 0.27 J ^ ' 

for e = 0.3. For other values of e in the range 0.2 to 0.35, the a„ oc ^ approximate law holds 
also, with the overall scale being e-dependent. 

Since the caustic rings are in or near the galactic plane, they cause bumps in the rotation 
curve at r ~ a„. Galactic rotation curves have bumps for many reasons unrelated to caustics. 
However the bumps caused by caustics have a special pattern. Eg. (14. II) . One may hope to 
find evidence for this pattern in a statistical analysis of many galactic rotation curves. 



In refs. [18|, Il9| 32 extended and well measured galactic rotation curves were selected 



under the criteria that each is an accurate tracer of the radial force law and that it extends 



far beyond the edge of the luminous disk. The data set of refs. [18|, ll9[ was analyzed as 



follows fl^. For each rotation curve, the radial variable was rescaled according to 

'220 km/s' 



"Vrot 



(4.2) 



where fj-ot is the rotation velocity read off from that particular curve. The data points for 
f < 10 kpc were deleted to remove the effect of the luminous disk. The 32 rotation curves 
were then co-added to make a combined rotation curve. 

There are two peaks in the combined rotation curve, one near 40 and one near 20 kpc, 
with significance of 3. Oct and 2.6(7, respectively. No explanation has been given for the 
occurence of the two peaks other than that they are the effect of the n = 1 and n = 2 
caustic rings of dark matter. The result suggests not only the existence of caustic rings of 
dark matter distributed according to Eq. (14. ip . but also that the jmax distribution of the 32 
galaxies in the data set is peaked near jmax.oid = 0.27 for e = 0.3. 

That the jmax distribution is peaked near jmax.ow = 0.27 was not a prediction of the 
self-similar infall model but it is certainly an interesting outcome of the analysis. It was 
mentioned in the previous subsection that the average magnitude j of the dimensionless 
angular momentum of the Milky Way had been determined to be jow ~ 0.2 from an estimate 
of the effective core radius of its halo. For a given j distribution on the turnaround sphere, 
jmax and i are related. The simplest j distribution is that of a rigidly rotating sphere 
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b(^) = imaxSin6'], in which case jmax = ^j- If we adopt this model, the estimate of the 
effective core radius of the Milky Way implies that its jmax,oid ~ ^ 0.2 = 0.255 . This is close 
to the peak in the jmax distribution of the 32 external galaxies. 

Finally, we mention that the rotation curve of NGC3198, which may be the best measured 
in the set [l^, 19l| . by itself shows three faint bumps at radii consistent with Eq. (14.11) and 

Jmax,old = 0.28 M- 



D. Milky Way rotation curve 

The effect of a caustic ring in the plane of a galaxy upon its rotation curve was analyzed 



in ref. [Uj. A caustic ring of radius a and width p produces a rise in the rotation curve which 
starts with an upward kink at ri = a and ends with a downward kink at r2 = a + p. The two 
discontinuities are a direct consequence of the fact that the dark matter density associated 
with the caustic ring diverges at r = a and r = a + p. 

For technical reasons, the Milky Way rotation curve is measured more precisely at r < r©, 
where is our own galactocentric distance, than at r > r©. We assume the standard value 
Tq = 8.5 kp c. The most detailed inner Galactic rotation curve, that we are aware of, was 



derived 23| from the Massachusetts-Stony Brook North Galactic Plane CO survey [241]. It 
has a series of sharp rises between 3 and 8.5 kpc [13]. Strikingly, each rise starts with an 
upward kink and ends with a downward kink, as expected for rises caused by caustic rings 
of dark matter. Where each rise starts (r = ri) and ends (r = r2), the slope of the rotation 
curve changes abruptly, from one data point to the next. 

Eq. (14. ip predicts ten caustic rings between 3 and 8.5 kpc, assuming jmax.oid — 0.255, 
as was inferred above from the effective core radius of the Milky Way halo. Allowing for 
ambiguities in identifying rises, the number of rises in the rotation curve [23;] between 3 and 
8.5 kpc is in fact ten plus or minus one. The radii where the observed rises start (^i „) and 
end {r2^n) are listed in Table HVl under entries n = 5,6 ... 14. The self-similar infall model 
predictions for the caustic ring radii were fitted to the rises l3] . For a given value of e, this 



is a one parameter fit, jmax being the only free parameter. For e = 0.3, the best fit occurs 
for jmax.oid = 0.263. This is consistent with the estimate joid ~ 0.2 obtained from fitting 
the effective core radius of the Milky Way (see subsections IV.B and IV. C). The root mean 
square relative deviation between the fitted caustic ring radii and the radii ri„, where 
the rises in the rotation curve start with an upward kink, is 3.1%. 

As explained in Section III. A, jmax.oid = 0.263 in Einstein - de Sitter cosmology corre- 
sponds to jmax — 0.179 in Concordance Cosmology. We repeated the fit of the calculated 
a„ to the observed ri „ in Concordance Cosmology and obtained jmax = 0.186, which is the 
value adopted here. The slight shift in the fitted value of jmax compared to ref. 12] is due to 



the fact that the a„ are only approximately proportional to jmax, so that the overall increase 
in length scales by the factor 1.47 (see Section IIA) between the Einstein - de Sitter and 
Concordance cosmologies is only approximately compensated for by dividing jmax by 1.47. 

As already mentioned, the rotation curve of the Milky Way is much less well measured for 
r > Tq. It does however have a prominent rise between 12.7 and 13.7 kpc, where Eq. (l4.ip 
predicts the n = 3 caustic ring to lie. See, for example, the Milky Way rotation curve 



published in ref. [25 
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E. Triangular feature in IRAS map 



The gravitational fields of caustic rings of dark matter may leave imprints upon the dis- 
tribution of ordinary matter. Looking tangentially to a caustic ring from a vantage point 
in the plane of the ring, one may recognize the tricusp shape of the D_4 catastrophe. The 
IRAS map of the galactic disk in the direction of galactic coordinates (Z, b) = (80°, 0°) shows 



a triangular shape which is strikingly suggestive of the cross-section of a caustic ring [12 
The relevant IRAS maps are posted at http://www.phys.ufl.edu/~sikivie/triangle/ . The 
triangular shape is correctly oriented with respect to the galactic plane and the galactic cen- 
ter. To an extraordinary degree of accuracy it is an isosceles triangle with axis of symmetry 
parallel to the galactic plane, as is expected for a caustic ring whose transverse dimensions 
are small compared to its radius. Moreover its position is consistent with the position of 
the rise in the rotation curve, between 8.28 and 8.43 kpc, caused by the caustic ring of dark 
matter nearest to us. 



F. The Monoceros Ring of stars 

The model predicts that the second caustic ring of dark matter (n = 2) lies in the 
Galactic plane at radius 02 — 20 kpc. After this prediction was made a ring of stars. 



named the "Monoceros Ring", was discovered in the Galactic plane at r ~ 20 kpc [26|. It is 



shown in ref. [27j that the Monoceros Ring of stars is the plausible outcome of the presence 
of the second caustic ring of dark matter. In particular, it was shown that the adiabatic 
modification of disk star orbits by the gravitational field of the caustic ring causes an order 
100% enhancement of the density of disk stars at the location of the caustic ring. It was 
also shown that viscous forces drive the gas in the neighborhood of the caustic ring towards 
the caustic ring radius 02 which may therefore be a site of enhanced star formation. 



V. THE FLOW DENSITIES AND VELOCITIES ON EARTH 

Our study of the phase space structure of the Milky Way halo is motivated in large part 
by the ongoing searches for dark matter on Earth, using axion m and WIMP Hi detectors. 



The signal in such detectors depends on the velocity distribution of dark matter in the solar 
neighborhood. The caustic ring halo model predicts that most of the local dark matter is 
in discrete flows, and predicts the velocity vectors and densities of the flows. Table |V] lists 
the properties of the first forty flows at the Earth's location in the Galaxy. Earlier versions 
of this table appeared in refs. jsoj and [lH. The purpose of this section is to describe how 
the estimates of Table |V] were obtained and to comment on the uncertainties. 

As described earlier, the flows come in pairs labeled (n, ±) . The flow velocities are stated 
in a reference frame which is attached to the Galaxy but which is not rotating relative to 
the faraway universe. The flow speeds (column 2) follow from the kinetic energy acquired 
by the particles while they fall in the growing gravitational potential well of the Galaxy. We 
calculate the speeds at the location of the Sun, at Tq = 8.5 kpc from the Galactic center, in 
the spherically symmetric self-similar model. Thus, 

Vn±{rQ) = J ^n± ^(^n±) (5.1) 
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where R = 2.4 Mpc, t = 13.7 Gyr, and the are the solutions of: 

re = Rm ■ (5.2) 

There are two solutions for each n but the corresponding speeds are nearly equal. The 
difference between Vn+{rQ) and ^^.(r©) is ignored in Table IVl 

The velocity components are listed in columns 3 through 5 of Table IVl and the densities in 
columns 6 and 7. The previous sections provide two ways to estimate the flow velocities and 
densities. For those flows whose inner caustic rings are close to us, we have the description 
of Section IIIC. For those flows whose inner caustic rings are far from us, we have the 
description of Section HID, and specifically the approximate equations fl3.45p . In either 
case, we assume that the Sun is in the symmetry plane [z = 0) of the model. This leads to 
some simplifications. 

For z = Eqs. (13.161) . which describe the flow near a caustic ring, have two pairs of 
solutions, corresponding to a = and r] = 0. For a = we have: 

r = p = a+ ^u{r]-T]oY 



dp / r 

^ _ 1 dM 1 

The choice of sign in the equation for Vp is independent of the choice of sign in the equation 
for d±. Eqs. (15.31) describe the in and out flows which exist in the plane of the ring for r > a. 
For ?7 = we have 

r = p = 
dz 

^ dr] 
dp 
df] 



Vp 




rdVLdt 2b{a+p-r) ' ^^'"^"^ 

The choice of sign in the equation for Vz is independent of the choice of sign in the equation 
for Vp. Eqs. (15. 4p describe the up and down flows which exist in the plane of the ring for 
r < a + p. 

For z = Eqs. (13.451) become 



(^v-±zvdl — ^- j for r < a 



= (fv— ±pvyl — y—j for r > a 

d -^"^^ ^ ' ^) 

~ vdQdt\r^-a^\ • ^^-^^ 
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The (only) justification for Eqs. (13.450 . and lience for Eqs. (15.51) . is that they provide a 
smooth interpolation of the flows between the inner and outer caustics in the limit where 
the transverse dimensions, p and q, of the inner caustic are neglected. For flows of n near 
5, whose caustic rings are near us, Eqs. (15. 3p and (15. 4p provide a more accurate description 
but they involve the parameters b and s about which we have little information other than 
that b is of order v, and s of order a. Eqs. (I3.45P and (15. 5p merely assume b = v and s = a. 
One may verify that, for b = v and s = a, Eqs. (15.50 are consistent with Eqs. (15. 3p and (15. 4p 
when r is near a but p = 0. 

The velocity components listed in columns 3 through 5 of Table IVl were obtained using 
Eqs. (15. 5p . For those values of n (i.e. ra = 3, 5 through 14) for which there is evidence for 
caustic rings of dark matter in the form of sharp rises in the Galactic rotation curve, we 
used the caustic ring radii ri „ inferred from the positions of the rises. The ri „ are listed in 
Table HVl For the other values of n (i.e. n = 1, 2, 4, 15 through 20) we used the predicted 
values listed in Table IIIII 

Eqs. (15. 5p predict that v^^Ivq) = for n > 5 and v^^Itq) = for n < A. That ^"^(r©) = 
for n > 5 follows from the symmetry of the model. We have set the corresponding entries 
equal to zero in Table IVl On the other hand, that Vp^{r^) = for n < 4 is a consequence 
of the approximations that led to Eqs. (13.450 . in particular the approximation p = 0. We 
may compare the zero value for Vp in Eqs. (15. 5p for r < a with the expression for Vp in Eqs. 
(15. 4p in terms of quantities v and a, which we believe we know well, and p on which there is 
observational information. The widths Pn of those caustic rings for which there is evidence 
in the form of sharp rises in the Galactic rotation curve may be obtained from Table IV by 
setting Pn = T2.n — Ti^n- FoT this set, the ratio | varies from 0.015 to 0.10, with an average 
of 0.05. Eqs. (15.40 imply therefore that, at least for r near a, the up and down flows for 
r < a + p have \vp\ ~ 0.3 v, whereas Eqs. (15. 5p would predict a zero value. As a reminder of 
these uncertainties, the entries in Table HVl for v'p^{rQ) with n < 5 have been left blank. 

The flow densities listed in the last two columns of Table IVl were calculated using Eqs. (15. 5p 
for n < 4 and n > 10, and using Eqs. (15. 3p with b = v for 5 < < 9. Again we used the 
observed values ri^n of the caustic ring radii for n = 3, 5 through 14 and the predicted for 
n = 1,2 4 and 15 through 20. The values of were taken from the last column of Table 
IIIII Since the ratio | was set equal to one, whereas it is only known to be of order one, all 
density estimates are uncertain by at least a factor two. 

VI. SUMMARY 

In this paper our main goal has been to provide a complete and self-contained description 
of the caustic ring model of the Milky Way halo, so that anyone wishing to compare the 
model with observations may readily obtain its predictions. Aside from the assumption of 
cold coUisionless dark matter, the defining properties of the model are axial and reflection 
symmetry, self-similarity and net overall rotation. The model is a proposal for the full phase 
space distribution of the Milky Way halo, which is described as a set of discrete flows with 
stated densities and velocities everywhere. 

In Section II, we discussed the general conditions under which the time evolution of a 
self-gravitating cold coUisionless fluid is self-similar. We found that self-similarity does not 
require symmetry. We discussed three cases: no symmetry, spherical symmetry, and axial 
symmetry. We derived the equations of motion for the reduced variables in each instance, 
and the initial conditions consistent with self-similarity. 
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The model breaks spherical symmetry because the dark matter particles carry angular 
momentum. The angular momentum distribution determines the properties of the inner 
parts of the halo, and in particular of the inner caustics. However, on large scales, the model 
is nearly spherically symmetric. The spherically symmetric model was originally described 
by Fillmore and Goldreich 0], and by Bertschinger Q]. It depends on one dimensionless 
parameter e and two dimensionful parameters. The latter merely set the scale of the halo's 
extent in phase space. They are determined in our model by requiring the Galactic rotation 
velocity at the location of the Sun to be 220 km/s and the age of the Galaxy to be 13.7 Gyr. 
The parameter e is related to the slope of the power spectrum of density perturbations on 
galactic scales. This implies e ~ 0.3. The properties of the model fitted to observation do 
not depend sharply on e. We set e = 0.3 for definiteness. 

The properties of the outer caustic spheres, the flow speeds and the mass infall rates 
were calculated in the spherically symmetric model. The properties of the first twenty outer 
caustics are listed in Table [TTl The flow speeds at the position of the Sun are listed in Table 
rVl The flow speeds at the location of each flow's inner caustic are listed in Table IIIII The 
mass infall rates are given by Eq. (13.411) . Their values at the inner caustic radii are given in 
Table [ml 

The structure of the inner caustics depends on the angular momentum distribution of 
the infalling dark matter particles. The additional parameters introduced in our axially 
symmetric and self-similar halo model, relative to the spherically symmetric model, are 
the components ^'r(a), I'eia) and z/^(q;) of the rescaled initial velocity distribution on the 
turnaround sphere, Eq. fl2.36p . The radial component //^(a) only plays a minimal role 
[sl because it does not contribute to angular momentum. We assume that the angular 
momentum distribution is such that there is net overall rotation, i.e. v^p{a) dominates over 
ve{(y)- We make this assumption because in that case the inner caustics are tricusp rings. 
There is evidence for such rings in the Milky Way and other isolated spiral galaxies. The 
evidence is summarized in Section IV. The radii of the caustic rings depend on a single 
parameter jmax = T-'ipiS^)- The evidence is consistent with jmax = 0.186, the value adopted 
here. The caustic ring radii are listed in Table lllli 

The flow velocities and densities in the neighborhood of a caustic ring are described in 
Section IIIC in terms of 6 parameters which may be chosen to be: the caustic ring radius 
a„, the flow speed Vn at the caustic, the mass infall rate ^^|n there, the time scale ?7o,n over 
which the particles traverse the caustic, a parameter 6„ of order and a parameter s„ of 
order an- The model values of the first three parameters are given in Table UTTl for the first 
twenty inner caustics. The last three parameters (?7o,n, &n and s„) are related to relatively 
subtle details of the initial velocity distribution, as described in subsection III.C.4. They 
determine the transverse sizes p„ and g„ of the caustic rings (see Eq. (I3.19P ). and also enter 
the formula for the density Eq. (13.171) . 

The model properties are obtained by solving a couple of equations (e.g. Eq. (13.221) ) 
on a computer. However, quick estimates can be gotten by using the method of adiabatic 
invariants. The main results of that method are stated in Eqs. (13.461) and (13.471) . The 
method of adiabatic invariants usually agrees with the results of numerical integration to 
within 10% or 20%. It allows estimates of the model properties for e 7^ 0.3. 

Table |V] lists the predicted densities and velocity vectors of the first forty flows on Earth. 
Knowledge of the local dark matter velocity distribution is essential when interpreting a 
signal in a dark matter detector on Earth. A striking property of the model is the existence 
on Earth of a "Big Flow", entry 5+ in Table IVl whose existence should become readily 
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apparent if a signal is found. The model makes predictions for all approaches to the dark 
matter problem, including direct searches for axions 0, [m and WIMPs |3l|,j32|, strong 
and weak gravitational lensin g [33| , WIMP annihilation in the Galactic halo j34| and the 
cosmic rays produced thereby jSSi ]. 
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TABLE I: The dimensionless rotation velocity squared 
turnaround radius R and the total mass M for different va' 



e 




R (Mpc) 


M(Mo) 


0.1 


0.21 


6.0 


3.1 X 10^^ 


0.15 


0.51 


3.9 


8.4 X 10^3 


0.2 


0.81 


3.1 


4.3 X 10^3 


0.25 


1.08 


2.7 


2.8 X 10^3 


0.3 


1.35 


2.4 


2.0 X 10^3 


0.35 


1.62 


2.2 


1.5 X 10^3 


0.4 


1.91 


2.0 


1.2 X 10^^ 


0.45 


2.15 


1.9 


9.8 X 10^2 



at the Sun's position i/^(e, the 
ues of e. 



TABLE IL The radii Rn and fold coefficients An of the first 20 outer caustics. 



ft 






1 

± 




7 7 V 1 n^^ 

I.I A lU 


2 


260 


8.8 X 10^^ 


3 


187 


9.6 X 10-^ 


4 


147 


1.02 X 10-=^ 


5 


121 


1.07 X 10-^ 


6 


103 


1.12 X 10^3 


7 


89 


1.16 X 10-3 


8 


79 


1.20 X 10-3 


9 


71 


1.24 X 10^3 


10 


64 


1.27 X 10-3 


11 


59 


1.29 X 10-3 


12 


54 


1.32 X 10-3 


13 


50 


1.35 X 10-3 


14 


47 


1.37 X 10-3 


15 


44 


1.40 X 10-3 


16 


41 


1.42 X 10-3 


17 


39 


1.44 X 10-3 


18 


37 


1.46 X 10-3 


19 


35 


1.48 X 10-3 


20 


34 


1.50 X 10-3 
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TABLE III: The radii a^, flow speeds Vn and infall rates ,„ ■ |„ for the flrst 20 inner caustics. 
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TABLE IV: Galactocentric radii where observed rises in the Milky Way rotation curves start (ri,n) 
and end {r2,n)- 



n 


ri,„ (kpc) 


r2,n (kpc) 


3 


12.7 


13.7 


5 


8.28 


8.43 


6 


7.30 


7.42 


7 


6.24 


6.84 


8 


5.78 


6.01 


9 


4.91 


5.32 


10 


4.18 


4.43 


11 


3.89 


4.08 


12 


3.58 


3.75 


13 


3.38 


3.49 


14 


3.16 


3.25 
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TABLE V: Velocity vectors iP^^{rQ) and densities ^^(r©) of the first 40 flows at our location in 
the Milky Way, in galactic coordinates, ip is in the direction of galactic rotation, z is perpendicular 
to the galactic plane, and p is in the radial direction. The two flows with velocity vector have 
density or d^, for each n independently. The Vp^{rQ) entries for n < 4 are not predicted but 
are expected to be relatively small; see text. 
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log(r) 



FIG. 1: Qualitative description of the phase space distribution of dark matter particles in a halo at 
a fixed moment of time. The solid lines represent occupied phase space cells. The vertical dotted 
line gives the observer position. Each intersection of the solid and dotted lines corresponds to a 
dark matter flow at the observer position. 
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FIG. 2: The function A(^) in the spherically symmetric self-similar infall model with e = 0.3. 
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FIG. 3: The phase space distribution of halo dark matter particles at a fixed moment of time in 
the spherically symmetric self-similar infall model with e = 0.3. The solid lines represent occupied 
phase space cells. 
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FIG. 4: Rotational velocity squared curves for different values of e, and j = 0. 
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FIG. 5: Dark matter trajectories forming a caustic ring of dark matter, in p-z cross-section. 




FIG. 6: The envelope of the trajectories that are shown in Fig. [5l The density diverges at the 
envelope in the limit of zero velocity dispersion. We refer to the shape shown as the "tricusp". 
What is meant by the radius a and the transverse sizes p and q oi a caustic ring is indicated. The 
tricusp has a discrete symmetry involving a rescaling and rotation by 120° about the central point, 
indicated by a star. 
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